clear;clc;
global theta e12 eg ey;

Nblack = 607;
Nwhite = 3973;

rng(10010);
N_sim=100000;

wparam=[0.174027339165297;0.187646010085453;1.52703723856717;0.578390833256311;1.47397496404516;0.472254474029496;0.539039975907426;0.508222124622837;0.557092030919351;1.68356685222368;-0.456906445371444;1.53322800556343;0.520506916649071;0.554621783550490;0.495990534196635;0.550405334264551;0.503398932635956;0.149831038482197;1.21750922668894;0.774068381366255;-0.0943955231130527;-0.454412676623209;0.225964287920840;-7.47819116937451e-05;0.226994852644431;0.162646835240202];
bparam=[-0.752900651375512;-0.625236744850452;0.955164978663548;0.470520928662698;0.940898873680747;0.457615100857030;0.522441225000711;0.796501681223757;0.526383943012879;1.38074364109675;-0.825533240871413;1.01277501698395;0.503940752081312;0.230058029470066;0.270326488861762;0.435689411502657;0.141913310661526;-0.625025361512185;0.665429321866805;0.834583035216824;-0.257622993708532;-0.213896500811884;-0.532316091400443;-0.518906855782177;0.0476680272326578;0.309398571910553];


%% Distribution of \theta
  
theta = normrnd(0,1,N_sim,1);
w.theta = theta*wparam(8);
b.theta = theta*bparam(8);
w.eg = normrnd(0,1,N_sim,1)*wparam(20);
b.eg = normrnd(0,1,N_sim,1)*bparam(20);
ey = normrnd(0,1,N_sim,1);

% cm = param(1);
% cr = param(2);
% phir = param(3);
% c1 = param(4);
% phi1 = param(5);
% sigem = param(6);
% siger = param(7);
% sigtheta = param(8);
% c2 = param(9);
% phi2 = param(10);
% cy = param(11);
% phim = param(12);
% gamma1 = param(13);
% gamma2 = param(14);
% betay = param(15);
% beta1 = param(16);
% beta2 = param(17);
% cg = param(18);
% phig = param(19);
% sigeg = param(20);
% c1d = param(21);
% phi1d = param(22);
% c2d = param(23);
% phi2d = param(24);
% beta1d = param(25);
% beta2d = param(26);


%% DIFFERENT CASES

%% Both White Teachers
param = bparam; theta = b.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
a=sim_74(param,1);
[f_bo,x_bo]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% Whites, SR
param = wparam; theta = w.theta; eg=w.eg;
a=sim_74(param,0);
[f_w,x_w]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% due to cy, theta and GPA
param = bparam; theta = w.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
param(11)=wparam(11);
param(4)=bparam(4)-bparam(11)+wparam(11);
param(9)=bparam(9)-bparam(11)+wparam(11);
param(18:20)=wparam(18:20);
a=sim_74(param,1);
[f_bwcytg,x_bwcytg]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% Outcome due to teacher expectations
param = bparam; theta = w.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
param(4:5)=wparam(4:5);
param(16)=wparam(16);
param(9:10)=wparam(9:10);
param(4)=wparam(4)-wparam(11)+bparam(11);
param(9)=wparam(9)-wparam(11)+bparam(11);
param(17)=wparam(17);
% cy param 
param(11)=wparam(11);
param(4)=bparam(4)-bparam(11)+wparam(11);
param(9)=bparam(9)-bparam(11)+wparam(11);
param(18:20)=wparam(18:20);
a=sim_74(param,0);
[f_bwte,x_bwte]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);


% Give BLACKS the WHITE PRODUCTION FUNCTION VALUES
param = bparam; theta = b.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
param(13:15)=wparam(13:15);
a=sim_74(param,1);
[f_bwy,x_bwy]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);



figure
subplot(2,2,1)
plot(x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks','Whites','location','southeast');

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);

subplot(2,2,2)
plot(x_bwcytg,f_bwcytg,'k.-',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);


subplot(2,2,3)
plot(x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G, T_j','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);


subplot(2,2,4)
plot(x_bwcytg,f_bwcytg,'k.-',x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta and G','Blacks,  White \theta, G, T_j.','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);

saveas(gcf,'white_teacher.jpg')





figure
plot(x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks','Whites','location','southeast');
xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);
saveas(gcf,'cf_1.jpg')

figure
plot(x_bwcytg,f_bwcytg,'k.-',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);
saveas(gcf,'cf_2.jpg')

figure

plot(x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G, T_j','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);
saveas(gcf,'cf_3.jpg')

figure



plot(x_bwcytg,f_bwcytg,'k.-',x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta and G','Blacks,  White \theta, G, T_j.','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);

saveas(gcf,'cf_4.jpg')
clear;
clc;

global theta e12 eg ey N_sim prdiff1b prdiff2b lowcutoff1 lowcutoff2 highcutoff1 highcutoff2 bparam wparam


%% Proportion of Math, English, Black, White
Ntotal = 3973+607;
Peb = 123/(84+123);
Pem = 100/164;
Nw = 3973; Nb = 607;
% init number of teachers
weor = 3775/3973; 
wmor = 3751/3973;

Nblack=607;
Nwhite = 3973;

%% 

rng(10010);

% theta = normrnd(0,1,Ntotal,1);

N_sim=100;


p=[-0.751980808619154,0.173675294186044;-0.624353919823129,0.187302821705864;0.968435023771807,1.55736791123648;0.477775564535539,0.556112164935507;0.775477992882744,1.32329370850257;0.447649758272938,0.463277982730318;0.520451846915934,0.536131631955990;0.786981589510950,0.499936033862615;0.573916339321584,0.541325133326204;1.15474203153640,1.55555645879034;-0.842336883069695,-0.472487903591744;1.03224424529235,1.57044922104757;0.129994542216687,0.546090110787827;0.654787908867107,0.574066902485244;0.280712856619634,0.499391157286847;0.485218918739055,0.579815459040526;0.297959269179098,0.524645829053341;-0.624397629713842,0.149556794871275;0.670177790543952,1.23221778368632;0.836160100864257,0.776509542464759;-0.256732070046633,-0.0405246290321471;-0.0981763407668954,-0.326246938973390;-0.583755273711407,0.204970642584455;-0.315411727796991,-0.0444559510729305;0.484129673862288,-0.102763672880495;-0.491223705626178,-0.0834210240868046;0.00764346578425141,0.186734752917938;0.123665148262387,0.209260916183449;0.442708378996374,0.430558964560301];

% p=[theta_bc, theta_wc];
pa=p;
pa(1:2,:)=p(13:14,:);pa(22:23,:)=p(25:26,:);
pa(3,:)=p(11,:);pa(5,:)=p(15,:); 
pa(4,:)=abs(p(8,:)); pa(21,:)=p(29,:);
pa(6,:)=p(4,:); pa(7,:)=p(9,:); pa(24,:)=p(21,:); pa(25,:)=p(23,:);
pa(10,:)=p(16,:); pa(11,:)=p(17,:);pa(28,:)=p(27,:);pa(29,:)=p(28,:);
pa(12,:)=p(1,:);pa(15,:)=p(12,:);pa(18,:)=p(6,:);
pa(13,:)=p(2,:);pa(16,:)=p(3,:);pa(19,:)=p(7,:);
pa(14,:)=p(18,:);pa(17,:)=p(19,:);pa(20,:)=p(20,:);
pa(8,:)=p(5,:); pa(9,:)=p(10,:);pa(27,:)=p(24,:);pa(26,:)=p(22,:);


p=pa;




Nsim = 100000;
rng(10010);
ey = normrnd(0,1,Nsim,1);

p(22:23,:)=0;
p(1,:)=[0.5118,0.5213];
p(2,:)=[0.2368,0.5562];
p(3,:)=[-0.8284,-0.4587];
p(5,:)=[0.2703,0.4985];



wparam = p(:,2);
bparam = p(:,1);

% 
% gamma1=param(1);gamma2=param(2);gamma12 = param(22);gamma22 = param(23);
% cy = param(3);betay = param(5);
% sigtheta = param(4);sig12 = param(21);
% c1 = param(6);c2 = param(7);c1d = param(24);c2d = param(25);
% phi1 = param(8);phi2 = param(9);phi1d = param(26);phi2d = param(27);
% beta1 = param(10);beta2 = param(11);beta1d = param(28);beta2d = param(29);
% cm = param(12);phim = param(15);sigm = param(18);
% cr = param(13);phir = param(16);sigr = param(19);
% cg = param(14);phig = param(17);sigg = param(20);

theta = linspace(-2*bparam(4),2*bparam(4),50);
theta=theta';

f=normpdf(theta,0,bparam(4));
figure
plot(theta,f,'k')
title('Distribution of Black \theta')
xlabel('\theta')
ylabel('pdf')
saveas(gcf,'fig0.jpg')










G=bparam(14)+bparam(17)*theta;
% Black Teacher, partial effect of T1, T2 at mean
PrY1E1Mm = normcdf(bparam(3)+(bparam(1))*(1-normcdf(bparam(3)+bparam(5)*G))...
    +(bparam(2))*(0.5300-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
PrY1E0Mm = normcdf(bparam(3)+(bparam(1))*(0-normcdf(bparam(3)+G))...
    +(bparam(2))*(0.5300-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
figure
plot(theta,PrY1E1Mm,'k.-',theta,PrY1E0Mm,'k')
ylim([0.05,0.45]);
legend('T1=1','T1=0')
title('Effect of Black English Teacher Expectation on Black Students')
xlabel('\theta')
ylabel('Probability')
saveas(gcf,'fig1.jpg')

% Black Teacher, partial effect of T2, T1 at mean
PrY1EmM1 = normcdf(bparam(3)+(bparam(1))*(0.4715-normcdf(bparam(3)+bparam(5)*G))...
    +(bparam(2))*(1-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
PrY1EmM0 = normcdf(bparam(3)+(bparam(1))*(0.4715-normcdf(bparam(3)+G))...
    +(bparam(2))*(0-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
figure
plot(theta,PrY1EmM1,'k.-',theta,PrY1EmM0,'k')
ylim([0.05,0.45]);
legend('T2=1','T2=0')
title('Effect of Black Math Teacher Expectation on Black Students')
xlabel('\theta')
ylabel('Probability')
saveas(gcf,'fig2.jpg')


% White Teacher, partial effect of T1, T2 at mean
PrY1E1Mmw = normcdf(bparam(3)+(bparam(1)+bparam(22))*(1-normcdf(bparam(3)+bparam(5)*G))...
    +(bparam(2)+bparam(23))*(0.5300-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
PrY1E0Mmw = normcdf(bparam(3)+(bparam(1)+bparam(22))*(0-normcdf(bparam(3)+G))...
    +(bparam(2)+bparam(23))*(0.5300-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
figure
plot(theta,PrY1E1Mmw,'k.-',theta,PrY1E0Mmw,'k')
ylim([0.05,0.45]);
legend('T1=1','T1=0')
title('Effect of White English Teacher Expectation on Black Students')
xlabel('\theta')
ylabel('Probability')
saveas(gcf,'fig3.jpg')

% White Teacher, partial effect of T2, T1 at mean
PrY1EmM1w = normcdf(bparam(3)+(bparam(1)+bparam(22))*(0.4715-normcdf(bparam(3)+bparam(5)*G))...
    +(bparam(2)+bparam(23))*(1-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
PrY1EmM0w = normcdf(bparam(3)+(bparam(1)+bparam(22))*(0.4715-normcdf(bparam(3)+G))...
    +(bparam(2)+bparam(23))*(0-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
figure
plot(theta,PrY1EmM1w,'k.-',theta,PrY1EmM0w,'k')
ylim([0.05,0.45]);
title('Effect of White Math Teacher Expectation on Black Students')
xlabel('\theta')
ylabel('Probability')
saveas(gcf,'fig4.jpg')
legend('T2=1','T2=0')

% Black Teacher, T1,T2 1 to 0
PrY1E1Mm = normcdf(bparam(3)+(bparam(1))*(1-normcdf(bparam(3)+bparam(5)*G))...
    +(bparam(2))*(1-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
PrY1E0Mm = normcdf(bparam(3)+(bparam(1))*(0-normcdf(bparam(3)+G))...
    +(bparam(2))*(0-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
figure
plot(theta,PrY1E1Mm,'k.-',theta,PrY1E0Mm,'k')
ylim([0.05,0.45]);
title('Effect of Black Teacher Expectations on Black Students')
xlabel('\theta')
ylabel('Probability')
saveas(gcf,'fig5.jpg')
legend('T1,T2=1','T1,T2=0')


% White Teacher, T1,T2 1 to 0
PrY1E1Mmw = normcdf(bparam(3)+(bparam(1)+bparam(22))*(1-normcdf(bparam(3)+bparam(5)*G))...
    +(bparam(2)+bparam(23))*(1-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
PrY1E0Mmw = normcdf(bparam(3)+(bparam(1)+bparam(22))*(0-normcdf(bparam(3)+G))...
    +(bparam(2)+bparam(23))*(0-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
figure
plot(theta,PrY1E1Mmw,'k.-',theta,PrY1E0Mmw,'k')
ylim([0.05,0.45]);
title('Effect of White Teacher Expectations on Black Students')
xlabel('\theta')
ylabel('Probability')
saveas(gcf,'fig6.jpg')
legend('T1,T2=1','T1,T2=0')

% White Eng, Black Math 1 to 0
PrY1E1Mmw = normcdf(bparam(3)+(bparam(1)+bparam(22))*(1-normcdf(bparam(3)+bparam(5)*G))...
    +(bparam(2))*(1-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
PrY1E0Mmw = normcdf(bparam(3)+(bparam(1)+bparam(22))*(0-normcdf(bparam(3)+G))...
    +(bparam(2))*(0-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
figure
plot(theta,PrY1E1Mmw,'k.-',theta,PrY1E0Mmw,'k')
ylim([0.05,0.70]);
title('Effect of Black English-White Math Teacher Expectations on Black Students')
xlabel('\theta')
ylabel('Probability')
saveas(gcf,'fig7.jpg')
legend('T1,T2=1','T1,T2=0')























% PrY1E1M1 = normcdf(bparam(3)+(bparam(1))*(1-normcdf(bparam(3)+bparam(5)*G))...
%     +(bparam(2))*(0-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
% PrY1E0M1 = normcdf(bparam(3)+(bparam(1))*(0-normcdf(bparam(3)+G))...
%     +(bparam(2))*(0-normcdf(bparam(3)+bparam(5)*G))+bparam(5)*G);
% 
% plot(theta,PrY1E1M1,theta,PrY1E0M1)







%%
c1tw = wparam(6)-wparam(3);
c2tw = wparam(7)-wparam(3);

Prt1w = normcdf(bparam(3)+c1tw+wparam(8)*theta+wparam(10)*G);
Prt1b = normcdf(bparam(6)+bparam(24)+(bparam(8)+bparam(26))*theta+(bparam(10)+bparam(28))*G);
% Prt1w = mean(Prt1w); Prt1b=mean(Prt1b);

figure
plot(theta,Prt1w,'k',theta,Prt1b,'k.-')
xlim([-2*bparam(4),2*bparam(4)]);
xlabel('\theta')
ylabel('Pr(T=1| \theta)')
legend('White Teacher, Debiased','White Teacher, Biased','location','southeast')
title('Pr(T=1| \theta, Black Student), White English Teacher')
saveas(gcf,'figure1.jpg')



prdiff1 = Prt1w-Prt1b;
Prt2w = normcdf(bparam(3)+c2tw+wparam(9)*theta+wparam(11)*G);
Prt2b = normcdf(bparam(7)+bparam(25)+(bparam(9)+bparam(27))*theta+(bparam(11)+bparam(29))*G);
% Prt2w = mean(Prt2w); Prt2b=mean(Prt2b);
prdiff2 = Prt2w-Prt2b;
figure
plot(theta,Prt2w,'k',theta,Prt2b,'k.-')
xlim([-2*bparam(4),2*bparam(4)]);
xlabel('\theta')
ylabel('Pr(T=1| \theta)')
legend('White Teacher, Debiased','White Teacher, Biased','location','southeast')
title('Pr(T=1| \theta, Black Student), White Math Teacher')
saveas(gcf,'figure2.jpg')





Prt1w = normcdf(bparam(6)+(bparam(8))*theta+(bparam(10))*G);
Prt1b = normcdf(bparam(6)+bparam(24)+(bparam(8)+bparam(26))*theta+(bparam(10)+bparam(28))*G);
% Prt1w = mean(Prt1w); Prt1b=mean(Prt1b);

figure
plot(theta,Prt1w,'k',theta,Prt1b,'k.-')
xlim([-2*bparam(4),2*bparam(4)]);
xlabel('\theta')
ylabel('Pr(T=1| \theta)')
legend('Black Teacher','White Teacher, Biased','location','southeast')
title('Pr(T=1| \theta, Black Student), English Teacher')
saveas(gcf,'figure3.jpg')



prdiff1 = Prt1w-Prt1b;
Prt2w = normcdf(bparam(7)+(bparam(9))*theta+(bparam(11))*G);
Prt2b = normcdf(bparam(7)+bparam(25)+(bparam(9)+bparam(27))*theta+(bparam(11)+bparam(29))*G);
% Prt2w = mean(Prt2w); Prt2b=mean(Prt2b);
prdiff2 = Prt2w-Prt2b;
figure
plot(theta,Prt2w,'k',theta,Prt2b,'k.-')
xlim([-2*bparam(4),2*bparam(4)]);
xlabel('\theta')
ylabel('Pr(T=1| \theta)')
legend('Black Teacher','White Teacher, Biased','location','southeast')
title('Pr(T=1| \theta, Black Student), Math Teacher')
saveas(gcf,'figure4.jpg')





%% Bias, All Sample

N_sim=100;

theta = normrnd(0,wparam(4),Nb,N_sim);
% theta = theta(theta>norminv(0.6));
% theta = unifrnd(0.6,1,Nb,N_sim);
% theta = norminv(theta,0,bparam(4));


Nb=Nblack;Nw=Nwhite;
% Nsim = 100000;

% ey = normrnd(0,1,Ntotal,N_sim);
b.ey = normrnd(0,1,Nb,N_sim);
w.ey = normrnd(0,1,Nw,N_sim);

b.sigtheta=bparam(4); w.sigtheta=wparam(4); b.sig12=bparam(21); w.sig12=wparam(21); b.sigg=bparam(20); w.sigg=wparam(20);

% % thetao = normrnd(0,1,Nsim,1);
% btheta = repmat(normrnd(0,1,Nb,1)*b.sigtheta,1,N_sim);
% wtheta = repmat(normrnd(0,1,Nw,1)*w.sigtheta,1,N_sim);
btheta = theta;

% 
% p1bb = 84/371;
% p2bb = 64/371;
% p1wb = 123/371;
% p2wb = 100/371;

b.e12 = mvnrnd([0,0],[1,b.sig12;b.sig12,1],Nb*N_sim);
b.e12a = b.e12(:,1);
b.e22a = b.e12(:,2);
b.e1 = reshape(b.e12a,Nb,[]);
b.e2 = reshape(b.e22a,Nb,[]);
b.e12=[b.e1,b.e2];


% eg =normrnd(0,1,Nsim,1);
b.eg = normrnd(0,1,Nb,N_sim)*b.sigg;
w.e12 = mvnrnd([0,0],[1,w.sig12;w.sig12,1],Nw*N_sim);
w.e12a = w.e12(:,1);
w.e22a = w.e12(:,2);
w.e1 = reshape(w.e12a,Nw,[]);
w.e2 = reshape(w.e22a,Nw,[]);
w.e12=[w.e1,w.e2];


% w.eg = eg*w.sigg;
w.eg= normrnd(0,1,Nw,N_sim)*w.sigg;



qr=linspace(0,Nb,5);
jr=linspace(0,Nb,50);

g=bparam(14)+bparam(17)*btheta+b.eg;



e12=b.e12; eg=b.eg; theta = btheta;ey=b.ey;


T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
kr = linspace(0,1,50);
ir = linspace(0,1,5);
for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
    T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,kr(i));
    ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
%     ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,i);
%     ab3(i,:,j)= simor4(bparam,T1wor,T2wor,3,i);
    
end
end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
% eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end









T1wor = 1-binornd(1,127/607,[Nb,N_sim]);
for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
%     T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,kr(i));
%     ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
    ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,kr(i));
%     ab3(i,:,j)= simor4(bparam,T1wor,T2wor,3,i);
    
end
end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end


T1wor = 1-binornd(1,100/607,[Nb,N_sim]);
T2wor = 1-binornd(1,127/607,[Nb,N_sim]);
% for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
%     T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     aw(i,:) = simor4(bparam,T1wor,T2wor,0,kr(i));
%     ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
%     ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,kr(i));
    ab3(i,:)= simor4(bparam,T1wor,T2wor,3,kr(i));
    
end
% end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end
eng3=mean(ab3,2);

figure
plot(kr,eng1(:,1),'k',kr,eng1(:,3),'k--',kr,eng1(:,5),'k.-')
title('Debiasing English Teachers')
legend('All White Eng. Teacher','50% White Eng Teacher','No White Eng. Techer')
ylim([0.23,0.34])
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'eng_allsample_s.jpg')

figure
plot(kr,eng2(:,1),'k',kr,eng2(:,3),'k--',kr,eng2(:,5),'k.-')
title('Debiasing Math Teachers')
legend('All White Math Teacher','50% White Math Teacher','No White Math Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
ylim([0.23,0.34])
saveas(gcf,'math_allsample_s.jpg')

figure
plot(kr,eng3,'k')
title('Debiasing All Teachers')
% legend('All White Teacher')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
ylim([0.23,0.34])
saveas(gcf,'all_allsample_s.jpg')



figure
plot(kr,eng1(:,1),'k',kr,eng1(:,3),'k--',kr,eng1(:,5),'k.-')
title('Debiasing English Teachers')
legend('All White Eng. Teacher','50% White Eng Teacher','No White Eng. Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'eng_allsample.jpg')

figure
plot(kr,eng2(:,1),'k',kr,eng2(:,3),'k--',kr,eng2(:,5),'k.-')
title('Debiasing Math Teachers')
legend('All White Math Teacher','50% White Math Teacher','No White Math Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'math_allsample.jpg')

figure
plot(kr,eng3,'k')
title('Debiasing All Teachers')
% legend('All White Teacher')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'all_allsample.jpg')














%% BIAS, top 2 quintile
N_sim=100;

% theta = normrnd(0,wparam(4),Nwhite,N_sim);
% theta = theta(theta>norminv(0.6));
theta = unifrnd(0.6,1,Nb,N_sim);
theta = norminv(theta,0,bparam(4));


Nb=Nblack;Nw=Nwhite;
% Nsim = 100000;

% ey = normrnd(0,1,Ntotal,N_sim);
b.ey = normrnd(0,1,Nb,N_sim);
w.ey = normrnd(0,1,Nw,N_sim);

b.sigtheta=bparam(4); w.sigtheta=wparam(4); b.sig12=bparam(21); w.sig12=wparam(21); b.sigg=bparam(20); w.sigg=wparam(20);

% % thetao = normrnd(0,1,Nsim,1);
% btheta = repmat(normrnd(0,1,Nb,1)*b.sigtheta,1,N_sim);
% wtheta = repmat(normrnd(0,1,Nw,1)*w.sigtheta,1,N_sim);
btheta = theta;

% 
% p1bb = 84/371;
% p2bb = 64/371;
% p1wb = 123/371;
% p2wb = 100/371;

b.e12 = mvnrnd([0,0],[1,b.sig12;b.sig12,1],Nb*N_sim);
b.e12a = b.e12(:,1);
b.e22a = b.e12(:,2);
b.e1 = reshape(b.e12a,Nb,[]);
b.e2 = reshape(b.e22a,Nb,[]);
b.e12=[b.e1,b.e2];


% eg =normrnd(0,1,Nsim,1);
b.eg = normrnd(0,1,Nb,N_sim)*b.sigg;
w.e12 = mvnrnd([0,0],[1,w.sig12;w.sig12,1],Nw*N_sim);
w.e12a = w.e12(:,1);
w.e22a = w.e12(:,2);
w.e1 = reshape(w.e12a,Nw,[]);
w.e2 = reshape(w.e22a,Nw,[]);
w.e12=[w.e1,w.e2];


% w.eg = eg*w.sigg;
w.eg= normrnd(0,1,Nw,N_sim)*w.sigg;



qr=linspace(0,Nb,5);
jr=linspace(0,Nb,50);

g=bparam(14)+bparam(17)*btheta+b.eg;



e12=b.e12; eg=b.eg; theta = btheta;ey=b.ey;


T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
kr = linspace(0,1,50);
ir = linspace(0,1,5);
for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
    T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,kr(i));
    ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
%     ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,i);
%     ab3(i,:,j)= simor4(bparam,T1wor,T2wor,3,i);
    
end
end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end









T1wor = 1-binornd(1,127/607,[Nb,N_sim]);
for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
%     T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,kr(i));
%     ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
    ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,kr(i));
%     ab3(i,:,j)= simor4(bparam,T1wor,T2wor,3,i);
    
end
end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end


T1wor = 1-binornd(1,100/607,[Nb,N_sim]);
T2wor = 1-binornd(1,127/607,[Nb,N_sim]);
% for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
%     T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     aw(i,:) = simor4(bparam,T1wor,T2wor,0,kr(i));
%     ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
%     ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,kr(i));
    ab3(i,:)= simor4(bparam,T1wor,T2wor,3,kr(i));
    
end
% end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end
eng3=mean(ab3,2);


figure
plot(kr,eng1(:,1),'k',kr,eng1(:,3),'k--',kr,eng1(:,5),'k.-')
title('Debiasing English Teachers, Top 2 Quintile')
legend('All White Eng. Teacher','50% White Eng Teacher','No White Eng. Techer')
ylim([0.49,0.61])
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'eng_s.jpg')

figure
plot(kr,eng2(:,1),'k',kr,eng2(:,3),'k--',kr,eng2(:,5),'k.-')
title('Debiasing Math Teachers, Top 2 Quintile')
legend('All White Math Teacher','50% White Math Teacher','No White Math Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
ylim([0.49,0.61])
saveas(gcf,'math_s.jpg')

figure
plot(kr,eng3,'k')
title('Debiasing All Teachers, Top 2 Quintile')
% legend('All White Teacher')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
ylim([0.49,0.61])
saveas(gcf,'all_s.jpg')



figure
plot(kr,eng1(:,1),'k',kr,eng1(:,3),'k--',kr,eng1(:,5),'k.-')
title('Debiasing English Teachers, Top 2 Quintile')
legend('All White Eng. Teacher','50% White Eng Teacher','No White Eng. Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'eng.jpg')

figure
plot(kr,eng2(:,1),'k',kr,eng2(:,3),'k--',kr,eng2(:,5),'k.-')
title('Debiasing Math Teachers, Top 2 Quintile')
legend('All White Math Teacher','50% White Math Teacher','No White Math Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'math.jpg')

figure
plot(kr,eng3,'k')
title('Debiasing All Teachers, Top 2 Quintile')
% legend('All White Teacher')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'all.jpg')















%% BIAS, 2nd
N_sim=100;

% theta = normrnd(0,wparam(4),Nwhite,N_sim);
% theta = theta(theta>norminv(0.6));
theta = unifrnd(0.6,0.8,Nb,N_sim);
theta = norminv(theta,0,bparam(4));


Nb=Nblack;Nw=Nwhite;
% Nsim = 100000;

% ey = normrnd(0,1,Ntotal,N_sim);
b.ey = normrnd(0,1,Nb,N_sim);
w.ey = normrnd(0,1,Nw,N_sim);

b.sigtheta=bparam(4); w.sigtheta=wparam(4); b.sig12=bparam(21); w.sig12=wparam(21); b.sigg=bparam(20); w.sigg=wparam(20);

% % thetao = normrnd(0,1,Nsim,1);
% btheta = repmat(normrnd(0,1,Nb,1)*b.sigtheta,1,N_sim);
% wtheta = repmat(normrnd(0,1,Nw,1)*w.sigtheta,1,N_sim);
btheta = theta;

% 
% p1bb = 84/371;
% p2bb = 64/371;
% p1wb = 123/371;
% p2wb = 100/371;

b.e12 = mvnrnd([0,0],[1,b.sig12;b.sig12,1],Nb*N_sim);
b.e12a = b.e12(:,1);
b.e22a = b.e12(:,2);
b.e1 = reshape(b.e12a,Nb,[]);
b.e2 = reshape(b.e22a,Nb,[]);
b.e12=[b.e1,b.e2];


% eg =normrnd(0,1,Nsim,1);
b.eg = normrnd(0,1,Nb,N_sim)*b.sigg;
w.e12 = mvnrnd([0,0],[1,w.sig12;w.sig12,1],Nw*N_sim);
w.e12a = w.e12(:,1);
w.e22a = w.e12(:,2);
w.e1 = reshape(w.e12a,Nw,[]);
w.e2 = reshape(w.e22a,Nw,[]);
w.e12=[w.e1,w.e2];


% w.eg = eg*w.sigg;
w.eg= normrnd(0,1,Nw,N_sim)*w.sigg;



qr=linspace(0,Nb,5);
jr=linspace(0,Nb,50);

g=bparam(14)+bparam(17)*btheta+b.eg;



e12=b.e12; eg=b.eg; theta = btheta;ey=b.ey;


T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
kr = linspace(0,1,50);
ir = linspace(0,1,5);
for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
    T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,kr(i));
    ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
%     ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,i);
%     ab3(i,:,j)= simor4(bparam,T1wor,T2wor,3,i);
    
end
end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end









T1wor = 1-binornd(1,127/607,[Nb,N_sim]);
for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
%     T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
    aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,kr(i));
%     ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
    ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,kr(i));
%     ab3(i,:,j)= simor4(bparam,T1wor,T2wor,3,i);
    
end
end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end


T1wor = 1-binornd(1,100/607,[Nb,N_sim]);
T2wor = 1-binornd(1,127/607,[Nb,N_sim]);
% for j = 1:length(ir);
for i = 1:length(kr);
%     T1wor = 1-binornd(1,123/607,[Nb,N_sim]);
%     T2wor = 1-binornd(1,100/607,[Nb,N_sim]);
%     T1wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     T2wor = 1-binornd(1,j/length(ir),[Nb,N_sim]);
%     aw(i,:) = simor4(bparam,T1wor,T2wor,0,kr(i));
%     ab(i,:,j) = simor4(bparam,T1wor,T2wor,1,kr(i));
%     ab2(i,:,j)= simor4(bparam,T1wor,T2wor,2,kr(i));
    ab3(i,:)= simor4(bparam,T1wor,T2wor,3,kr(i));
    
end
% end



for i=1:5
eng(:,i) = mean(aw(:,:,i),2);
eng1(:,i)=mean(ab(:,:,i),2);
eng2(:,i) = mean(ab2(:,:,i),2);
% eng3(:,i)=mean(ab3(:,:,i),2);
end
eng3=mean(ab3,2);


figure
plot(kr,eng1(:,1),'k',kr,eng1(:,3),'k--',kr,eng1(:,5),'k.-')
title('Debiasing English Teachers, 2nd Quintile')
legend('All White Eng. Teacher','50% White Eng Teacher','No White Eng. Techer')
ylim([0.35,0.5])
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'eng_s2.jpg')

figure
plot(kr,eng2(:,1),'k',kr,eng2(:,3),'k--',kr,eng2(:,5),'k.-')
title('Debiasing Math Teachers, 2nd Quintile')
legend('All White Math Teacher','50% White Math Teacher','No White Math Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
ylim([0.35,0.5])
saveas(gcf,'math_s2.jpg')

figure
plot(kr,eng3,'k')
title('Debiasing All Teachers, 2nd Quintile')
% legend('All White Teacher')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
ylim([0.35,0.5])
saveas(gcf,'all_s2.jpg')



figure
plot(kr,eng1(:,1),'k',kr,eng1(:,3),'k--',kr,eng1(:,5),'k.-')
title('Debiasing English Teachers, 2nd Quintile')
legend('All White Eng. Teacher','50% White Eng Teacher','No White Eng. Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'eng2.jpg')

figure
plot(kr,eng2(:,1),'k',kr,eng2(:,3),'k--',kr,eng2(:,5),'k.-')
title('Debiasing Math Teachers, 2nd Quintile')
legend('All White Math Teacher','50% White Math Teacher','No White Math Techer')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'math2.jpg')

figure
plot(kr,eng3,'k')
title('Debiasing All Teachers, 2nd Quintile')
% legend('All White Teacher')
xlabel('Percentage debiased')
ylabel('Pr(Y=1)')
saveas(gcf,'all2.jpg')





















% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% % figure;
% % 
% % plot(kr,eng(:,1),kr,eng(:,2),kr,eng(:,3),kr,eng(:,4),kr,eng(:,5))
% % legend('none','25','50','75','100')
% % 
% % plot(kr,eng1(:,1),kr,eng1(:,2),kr,eng1(:,3),kr,eng1(:,4),kr,eng1(:,5))
% % legend('none','25','50','75','100')
% % 
% % 
% % plot(kr,eng2(:,1),kr,eng2(:,2),kr,eng2(:,3),kr,eng2(:,4),kr,eng2(:,5))
% % legend('none','25','50','75','100')
% % 
% % 
% % plot(kr,eng3(:,1),kr,eng3(:,2),kr,eng3(:,3),kr,eng3(:,4),kr,eng3(:,5))
% % legend('none','25','50','75','100')
% 
% 
% 
% 
% 
% 
% 
% 
% 
% % plot1=plot(jr2,eng2(:,1),'k.-',jr2,eng2(:,3),'k',jr2,eng2(:,5),'k--',jr2,eng(:,1),'g.-',jr2,eng(:,3),'g',jr2,eng(:,5),'g--');
% 
% 
% 
% 
% 
% 
% 
% for j = 1:length(qr);
% for i = 1:length(jr);
%     T1wor = 1-binornd(1,jr(i)/Nb,[Nb,N_sim]);
%     T2wor = 1-binornd(1,qr(j)/Nb,[Nb,N_sim]);
%     e12=b.e12; eg=b.eg; theta = btheta;ey=b.ey;
%         aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,0);
%         ab(i,:,j)=simor4(bparam,T1wor,T2wor,1,1);
% end
% end
% 
% for i=1:5
% eng(:,i) = mean(aw(:,:,i),2);
% eng2(:,i)=mean(ab(:,:,i),2);
% end
% jr2 = jr/Nb*100;
% 
% 
% figure;
% 
% plot1=plot(jr2,eng2(:,1),'k.-',jr2,eng2(:,3),'k',jr2,eng2(:,5),'k--',jr2,eng(:,1),'g.-',jr2,eng(:,3),'g',jr2,eng(:,5),'g--');
% % plot1=plot(jr2,eng2(:,1),'k.-',jr2,eng2(:,5),'k--',jr2,eng(:,1),'b.-',jr2,eng(:,5),'b--');
% 
% xlabel('Percent of Black Students with Black Teachers')
% ylabel('Graduation Probability')
% title('Change in Graduation Probability by number of Black English Teachers')
% % legend('No Black Math Teachers, Debiased','50% Black Math Teachers, Debiased','All Black Math Teachers, Debiased',...
% %     'No Black Math Teachers','50% Black Math Teachers','All Black Math Teachers','location','southoutside')
% legend_str={'No Black Math Teachers, Debiased','50% Black Math Teachers, Debiased','All Black Math Teachers, Debiased',...
%     'No Black Math Teachers','50% Black Math Teachers','All Black Math Teachers'};
% % columnlegend(2,legend_str,'Location','SouthEast')
% 
% gridLegend(plot1,2,legend_str,'location','north','Fontsize',7,'Box','on');
% 
% saveas(gcf,'engblackdebias1.jpg')
% 
% 
% 
% 
% 
% 
% 
% figure
% plot(jr2,eng(:,1),'k.-',jr2,eng(:,3),'k',jr2,eng(:,5),'k--');
% xlabel('Percentage of Black Students with Black Teachers')
% ylabel('Proportion Graduating from 4-Year College')
% title('Change in Graduation Probability by number of Black English Teachers')
% legend('No Black Math Teachers','50% Black Math Teachers','All Black Math Teachers','location','northeast')
% saveas(gcf,'engblack1.jpg')
% 
% figure
% plot(jr2,eng(:,1),'k.-',jr2,eng(:,3),'k',jr2,eng(:,5),'k--');
% xlabel('Percentage of Black Students with Black Teachers')
% ylabel('Proportion Graduating from 4-Year College')
% title('Change in Graduation Probability by number of Black English Teachers')
% ylim([0,0.5]);
% legend('No Black Math Teachers','50% Black Math Teachers','All Black Math Teachers','location','southeast')
% saveas(gcf,'engblack2.jpg')
% 
% 
% 
% for j = 1:length(qr);
% for i = 1:length(jr);
%     T1wor = 1-binornd(1,qr(j)/Nb,[Nb,N_sim]);
%     T2wor = 1-binornd(1,jr(i)/Nb,[Nb,N_sim]);
%     e12=b.e12; eg=b.eg; theta = btheta;ey=b.ey;
%         aw(i,:,j) = simor4(bparam,T1wor,T2wor,0,0);
%                 ab(i,:,j)=simor4(bparam,T1wor,T2wor,2,1);
% 
% end
% end
% 
% for i=1:5
% math(:,i) = mean(aw(:,:,i),2);
% math2(:,i)=mean(ab(:,:,i),2);
% end
% jr2 = jr/Nb*100;
% 
% 
% figure
% plot2=plot(jr2,math2(:,1),'k.-',jr2,math2(:,3),'k',jr2,math2(:,5),'k--',jr2,math(:,1),'g.-',jr2,math(:,3),'g',jr2,math(:,5),'g--');
% 
% xlabel('Percent of Black Students with Black Teachers')
% ylabel('Graduation Probability')
% title('Change in Graduation Probability by number of Black Math Teachers')
% % legend('No Black Eng Teachers, Debiased','50% Black Eng Teachers, Debiased','All Black Eng Teachers, Debiased',...
% %     'No Black Eng Teachers','50% Black Eng Teachers','All Black Eng Teachers','location','northeast')
% legend_str={'No Black Eng Teachers, Debiased','50% Black Eng Teachers, Debiased','All Black Eng Teachers, Debiased',...
%     'No Black Eng Teachers','50% Black Eng Teachers','All Black Eng Teachers'};
% % columnlegend(2,legend_str,'Location','SouthEast')
% 
% gridLegend(plot2,2,legend_str,'location','north','Fontsize',7,'Box','on');
% 
% 
% 
% saveas(gcf,'mathblackdebias1.jpg')
% 
% 
% 
% 
% 
% 
% figure
% plot(jr2,math(:,1),'k.-',jr2,math(:,3),'k',jr2,math(:,5),'k--');
% xlabel('Percentage of Black Students with Black Teachers')
% ylabel('Proportion Graduating from 4-Year College')
% title('Change in Graduation Probability by number of Black Math Teachers')
% legend('No Black English Teachers','50% Black English Teachers','All Black English Teachers','location','northeast')
% saveas(gcf,'mathblack1.jpg')
% 
% 
% figure
% plot(jr2,math(:,1),'k.-',jr2,math(:,3),'k',jr2,math(:,5),'k--');
% xlabel('Percentage of Black Students with Black Teachers')
% ylabel('Proportion Graduating from 4-Year College')
% title('Change in Graduation Probability by number of Black Math Teachers')
% ylim([0,0.5]);
% legend('No Black English Teachers','50% Black English Teachers','All Black English Teachers','location','southeast')
% saveas(gcf,'mathblack2.jpg')
% 
% 
% 
% 
% %%
% % qr=linspace(0,Nb,5);
% % jr=linspace(0,Nb,50);
% % for j = 1:length(qr);
% for i = 1:length(jr);
%     T1wor = 1-binornd(1,jr(i)/Nb,[Nb,N_sim]);
%     T2wor = 1-binornd(1,jr(i)/Nb,[Nb,N_sim]);
%     e12=b.e12; eg=b.eg; theta = btheta;ey=b.ey;
%         Qaw(i,:) = simor4(bparam,T1wor,T2wor,0,0);
%                 Qab(i,:)=simor4(bparam,T1wor,T2wor,3,1);
% 
% end
% % end
% 
% % for i=1:5
% Amath = mean(Qaw,2);
% Amath2=mean(Qab,2);
% % end
% jr2 = jr/Nb*100;
% 
% 
% figure
% plot(jr2,Amath,'k');
% xlabel('Percent of Black Students with Black Teachers')
% ylabel('Graduation Probability')
% title('Change in Graduation Probability by number of Black Teachers')
% saveas(gcf,'figure_1.jpg')
% 
% 
% 
% 
% 
% figure
% plot2=plot(jr2,Amath2,'k.-',jr2,Amath,'k--');
% xlabel('Percent of Black Students with Black Teachers')
% ylabel('Graduation Probability')
% title('Change in Graduation Probability by number of Black Teachers')
% legend('Debiased','Baseline')
% % legend('No Black Eng Teachers, Debiased','50% Black Eng Teachers, Debiased','All Black Eng Teachers, Debiased',...
% %     'No Black Eng Teachers','50% Black Eng Teachers','All Black Eng Teachers','location','northeast')
% % legend_str={'No Black Eng Teachers, Debiased','50% Black Eng Teachers, Debiased','All Black Eng Teachers, Debiased',...
% %     'No Black Eng Teachers','50% Black Eng Teachers','All Black Eng Teachers'};
% % columnlegend(2,legend_str,'Location','SouthEast')
% 
% % gridLegend(plot2,2,legend_str,'location','north','Fontsize',7,'Box','on');
% 
% 
% 
% saveas(gcf,'Amathblackdebias1.jpg')
% 
% 
% 
% 
% % 
% % 
% % figure
% % plot(jr2,math(:,1),'k.-',jr2,math(:,3),'k',jr2,math(:,5),'k--');
% % xlabel('Percentage of Black Students with Black Teachers')
% % ylabel('Proportion Graduating from 4-Year College')
% % title('Change in Graduation Probability by number of Black Math Teachers')
% % legend('No Black English Teachers','50% Black English Teachers','All Black English Teachers','location','northeast')
% % saveas(gcf,'mathblack1.jpg')
% % 
% % 
% % figure
% % plot(jr2,math(:,1),'k.-',jr2,math(:,3),'k',jr2,math(:,5),'k--');
% % xlabel('Percentage of Black Students with Black Teachers')
% % ylabel('Proportion Graduating from 4-Year College')
% % title('Change in Graduation Probability by number of Black Math Teachers')
% % ylim([0,0.5]);
% % legend('No Black English Teachers','50% Black English Teachers','All Black English Teachers','location','southeast')
% % saveas(gcf,'mathblack2.jpg')
% 
% 
% 
% 
% 
% % black = mean(ab,2);
% % white = mean(aw,2);
% 
% 
% 
% 
% % 
% % 
% % for i = 1:2000;
% %         T1bor= 1-binornd(1,p1b(i),[Nb,N_sim]);
% %         T2bor= 1-binornd(1,p2b(i),[Nb,N_sim]);
% %         T1wor= binornd(1,p1w(i),[Nw,N_sim]);
% %         T2wor= binornd(1,p2w(i),[Nw,N_sim]);            
% %         e12=b.e12; eg=b.eg; theta = btheta;ey=b.ey;
% %         ab(i,:) = simor2(bparam,T1bor,T2bor,0);
% %         e12=w.e12; eg=w.eg; theta = wtheta;ey=w.ey;
% %         aw(i,:) = simor2(wparam,T1wor,T2wor,0);
% %     i
% % end
% % 
% % black = mean(ab,2);
% % white = mean(aw,2);
% % 
% % 
% % 
% % 
% % 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% % gamma1=param(1);gamma2=param(2);gamma12 = param(22);gamma22 = param(23);
% % cy = param(3);betay = param(5);
% % sigtheta = param(4);sig12 = param(21);
% % c1 = param(6);c2 = param(7);c1d = param(24);c2d = param(25);
% % phi1 = param(8);phi2 = param(9);phi1d = param(26);phi2d = param(27);
% % beta1 = param(10);beta2 = param(11);beta1d = param(28);beta2d = param(29);
% % cm = param(12);phim = param(15);sigm = param(18);
% % cr = param(13);phir = param(16);sigr = param(19);
% % cg = param(14);phig = param(17);sigg = param(20);
% 
% % t1black_star = normcdf(bparam(6)+bparam(8)*theta+bparam(10)*G);
% % t2black_star = normcdf(bparam(7)+bparam(9)*theta+bparam(10)*G);
% % t1dbwhite_star = normcdf(wparam(6)+wparam(8)*theta+wparam(10)*G);
% % t2dbwhite_star = normcdf(wparam(7)+wparam(9)*theta+wparam(10)*G);
% % t1whtie_star = normcdf(bparam(6)+bparam(24)+(bparam(8)+bparam(26))*theta+(bparam(10)+bparam(28))*G);
% % t2whtie_star = normcdf(bparam(6)+bparam(24)+(bparam(8)+bparam(26))*theta+(bparam(10)+bparam(28))*G);
% 
% % 
% % figure
% % plot(theta+bparam(3),prdiff1,'k',theta+bparam(3),prdiff2,'k.-')
% % legend('English Teacher','Math Teacher')
% % ylim([-0.4,0.4]);
% % xlim([bparam(3)-2*bparam(4),bparam(3)+2*bparam(4)]);
% % hold on;
% % plot(xlim,[0,0],'k--');
% % title('\Delta Pr(T=1|\theta) for White and Black Students, White Teachers')
% % xlabel('c_y+\theta')
% % ylabel('\Delta in Probability')
% % 
% 
% 
% 
% 
% 
% 
% 
